First-principles approach to rotational-vibrational frequencies and infrared intensity 

for H 2 adsorbed in nanoporous materials 
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The absorption sites and the low-lying rotational and vibrational (RV) energy states for H2 ad- 
sorbed within a metal-organic framework are calculated via van der Waals density functional theory. 
The induced dipole due to bond stretching is found to be accurately given by a first-principles driven 
approximation using maximally-localized- Wannier-function analysis. The strengths and positions 
of lines in the complex spectra of RV transitions are in reasonable agreement with experiment, and 
in particular explain the experimentally mysteriously missing primary line for para hydrogen. 

PACS numbers: 68.43.Bc, 78.30.-j, 82.75.-z 



Gas adsorption into nanoporous materials is of great 
interest for both fundamental science and applications. 
Molecular H2 is challenging because it can vibrate, ro- 
tate, and translate quantum mechanically about its bind- 
ing site due to its small mass. The vibration-rotation 
(RV) excitations induced by infrared (IR) absorption 
thus provide rich information However, determin- 
ing the origin and strength of these lines is challeng- 
ing because large unit cells are encountered in typical 
nanoporous structures, and the dynamic dipole is dis- 
tributed over spatially remote parts of the structure. To 
determine the absorption intensity, a precisely tractable 
experimental quantity, one must not only calculate the 
dipole, but also evaluate the quantum mechanical matrix 
element. An effective approximation scheme for doing 
this has not hitherto been found. 

Here, we present such a scheme based on the combi- 
nation of a self-consistent van der Waals density func- 
tional (vdW-DF) approach with maximally-localized- 
Wannier-function (MLWF) analysis 0, 0| and apply it 
to H2 adsorption in a prototypical metal-organic frame- 
work, MOF-5 Such materials have been extensively 
explored for hydrogen storage [6|], gas separation, catal- 
ysis, and sensors [7j. We analyze the dynamical proper- 
ties of the adsorbed H2, finding results consistent with 
experiment. Importantly, we apply the MLWF analy- 
sis to calculate the induced dipole moment due to H2 
adsorption and bond stretching, decomposing the dipole 
into the contributions from both adsorbed dihydrogen 
and MOF. Monitoring the change in each Wannier cen- 
ter of the MOF structure upon H2 adsorption provides 
an intuitive picture by breaking the H2-sorbent interac- 
tion into individual components of the MOF structure, 
thus identifying the parts that directly interact with the 
dihydrogen. Such knowledge is important to optimize 
MOF structures for desired properties. In the present 
case, we use this information to calculate the dynamical 
dipole moment and its matrix element for H2 vibrational 
transitions and RV transitions. We find that the IR in- 
tensity of the purely vibrational mode for para hydrogen 



is only about 2.5% of that for ortho hydrogen at the pri- 
mary adsorption site, which agrees beautifully with the 
missing line in the experiment. [f|. A selection rule for 
RV transitions at the relevant site is also obtained and 
supported by the IR data. 

The H2 binding sites are efficiently determined by self- 
consistent vdW-DF calculations 2]. A series of total 
energy calculations for different bond lengths, orienta- 
tions, and center-of-mass positions respectively are per- 
formed keeping the MOF atoms fixed at experimental 
positions [9j- The resulting potential energy surfaces are 
then used in the corresponding radial and rigid rotor 
Schrodinger equations respectively to extract the vibra- 
tional, rotational and translational frequencies 13, 
Anharmonic effects are fully included. 

It has been shown that the sum of the Wannier- 
function centers is connected to the Berry phase theory of 
bulk polarization [3|. The dipole in the unit cell is given 
by u = eY, m Z m R m - eJ2 n ,s P in r n, where Z m and R m 
are the atomic number and position of the m th nuclei 
and r n is the center of the n th Wannier function. Im- 
portantly, it is trivial to decompose the total dipole into 
components in various parts of the structure [J], which 
goes beyond the Berry-phase method. Thereby, we may 
use the change of Wannier center upon adsorption as a 
qualitative measure for understanding the H2-MOF in- 
teraction and to determine the important parts of the 
MOF that directly interact with hydrogen. 

There are four types of adsorption sites in this struc- 
ture, as established experimentally [l2T ] and theoreti- 
cally with reasonable agreement. We start with the 
positions determined by neutron scattering [12j and re- 
lax the H 2 with the vdW-DF approach, thereby confirm- 
ing the positions of the four sites, named the cup, 03, 
02, and benzene sites [l2j. Fig. Q] shows the position 
of the cup site and a portion of the MOF-5 structure 
where there exists 3-fold rotation symmetry among the 
three benzene ring branches. The distance between the 
H2 center-of-mass position and the oxygen atom passing 
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FIG. 1: The primary adsorption site and the change of Wan- 
nier centers (blue balls) due to H2 adsorption compared to 
bare MOF and free H2. Vector lengths are enlarged by 1200. 
The structure shown is just a fraction of the unit cell. 



TABLE I: Theory vs experiment [H for the stretch frequency 
shift of the adsorbed H2 relative to free H2 . See text for zero 
point energies (not included in the binding energies E_b here). 



site 


Theory 


Expr. 


Calculated E_b 




(cm" 1 ) 


(cm" 1 ) 


(kJ/mol) 


cup 


-23 


-27.5 


-11.1 


02 


-22 


-19.0 


-7.9 


03 


-13 


-17 


-7.8 


benzene 


-15 




-5.4 



through the rotation axis is about 4. 2 A which is some- 
what larger than the measured value of 3. 8 A [12j due to 



a known vdW-DF overestimation of bond lengths (14 1. 
Fig. [T] also shows the shift of the Wannier centers upon 
H2 adsorption with respect to the bare MOF and the 
free H2. See Supplemental Material for other sites 21 1. 
These figures show that the Wannier centers associated 
with the 7r bonds in the benzene ring change significantly 
upon H2 adsorption for all four adsorption sites, show- 
ing a clear and intuitive picture of the MOF components 
that interact directly with the adsorbed H2. 

Tabic U shows both the theoretical and experimental 
stretch frequency shifts of the adsorbed H2 with respect 
to the corresponding free H2 value. The agreement is 
good, within a reasonable error. Importantly, the origin 
of the IR peaks at —19 and —17 cm- 1 (not understood 
in the experiment [H) is now unraveled with the aid of 
our calculational results. Note that the calculated bind- 
ing energies at 02 and 03 sites are very close. However, 
vdW-DF t ypic ally overestimates intermediate-range in- 
teractions [15| , Since 03 has three benzene neighbors 
while 02 has two, the overestimation for the 03 site is 
expected to be larger than that for the 02 site. As a 
result, the 02 site is probably more favorable energeti- 
cally and should get populated more than the 03 site. 
This is consistent with the measurements where the —17 
cm- 1 peak is quite weak and appears only as a shoulder 
to the main line at —19 cm -1 . We therefore assign the 
— 19 cm -1 peak to the 02 site whereas —17 cm -1 to 03. 
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FIG. 2: left: Orientational dependence of the binding energy 
(kJ/mol) at the cup site; right: Ground state \tp r ot\ 2 ■ The 
distance from the origin is the probability for that orientation. 



The IR spectra also show some RV lines where both 
vibrational and rotational states change during a single 
transition. Usually inelastic neutron scattering is em- 
ployed to study the H2 rotational states and has been 
already applied to H2 in MOF-5 [5j, |6j. However, the low 
energy resolution limits detailed analysis. We therefore 
consider the RV measured with IR [fj for comparison 
with our rotational calculations. The left panel of Fig. [5] 
is the angular potential energy surface at the cup site. 
The coordinate system is chosen so that the origin is at 
the cup site and the Z axis is the 3-fold rotation axis (see 
Fig. S2 in Supplemental Material [7]). Fig. [2] shows that 
H2 tends to lie in the XY plane and to be perpendic- 
ular to the rotation axis (Z). The energies for in-plane 
orientations are almost uniform. Therefore, the rotation 
is essentially two dimensional, as shown by the flattened 
ground-state angular wave function in the right panel 
of Fig. [21 Combining the stretch frequency and the ro- 
tational energies (see Supplemental Material), we obtain 
the RV frequencies. The results for the cup site are shown 
as S transitions in Table HT1 where the frequency shifts are 
listed relative to the corresponding free H2 values (see 
Supplemental Material for other sites). The magnitude 
of the shifts is consistent between theory and experiment, 
particularly for the leading peaks in each category that 
are most intense. 

We also calculated the translational frequencies at the 
cup site associated with the motion of the whole H2 
against the adsorption site. The three translational fre- 
quencies, at 95, 108 and 133 cm -1 respectively, are con- 
sistent with the value of 84 cm -1 extracted from IR 
spectra [f|. They are also similar to that observed for 
H2 in Ceo (HO cm -1 ) 17]. The determination of the 
rotational and translational states gives the correspond- 
ing zero point energies of ~0.5 and 2 kJ/mol for H2 at 
the cup site. The binding energy after corrections is 
therefore about 8.5 kJ/mol and somewhat larger than 
the measured adsorption enthalpy of ^5 kJ/mol (l8l.[i"9j]. 
This overestimation by vdW-DF, also found in other 
MOF materials [lOj, is attributed to overestimation of 
the intermediate-range interactions [lij . 

The measured IR spectra for the cup adsorption site 
shows a strong pure vibrational peak due to the ortho- 
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H2, while the corresponding para line is not observed. 
Since the orientational energy map only shows a small 
rotational barrier, the missing para-H 2 line cannot be 
explained by the assumption of a frozen H 2 orientation. 
Moreover, the local structure around this site has Cs v 
symmetry. The rotational state of the para H2 has the 
same symmetry as Z and transforms as Al. Therefore 
the transitions between two Al states should be IR ac- 
tive, even though the X and Y components of the dipole 
give a vanishing contribution by symmetry. 

To understand the unexpected missing para-H2 line 
and to calculate the line weights in the more complex RV 
spectra, we evaluate the transition dipole integral explic- 
itly. Assuming the electronic state remains in the ground 
state and the RV wave function is separable, one has 
I a = (^lih^loMctWvib^lot)! where u is the dipole mo- 
ment and a = X, Y, or Z; the translational motion asso- 
ciated with H 2 center-of-mass is not included. The dipole 
is a function of the H2 internuclear distance, R, and the 
bond orientation is defined by (9, 0). It can be expanded 
as u a (R, 9, 0) — u a (Ro, 9, 0) + u' a (R, 9, <f>)\ Ro AR where 
Rq is the equilibrium bond length, and u' a is the deriva- 
tive of u a with respect to R. Since the vibrational wave 
functions depend only on the inter-nuclear distance, the 
integral of the first term vanishes for transitions between 
different vibrational states due to orthogonality. We 
find / Q = (rtjAR\r v . lb ){^LK(Ro,8^)\r rot }, where 
\tpr 0t ) = \jirrii) with j even (odd) for para (ortho) H2, 
and similarly for \ipl ot ). The radial integral is a constant 
for both ortho and para H2 and therefore unnecessary 
for understanding the missing line of para-H 2 . The an- 
gular integral determines the relative intensity between 
them. We now need to evaluate this integral, for which 
u' a (Ro,9,(f>) remains to be calculated. 

To perform ab initio calculations for u' a for every (#,0) 
is computationally expensive and impractical for this sys- 
tem. A possible approach is to compute the dipole from 
first principles for a few H2 orientations and derive from 
them the dipole of all the other orientations. This be- 
comes feasible if one can write the dipole as 

u a = J2 C i,<*FiA0,<t>) (!) 

i 

where F are some known functions and the summation 
needs to be run over only a few terms. This approach is 
appropriate if one realizes that H2 and MOF are weakly 
interacting and the dipole induced on each other can be 
well described within a classical picture. First, MOF 
atoms produce an electric field (E) which induces a dipole 
on H 2 . At the cup site, the field is along Z due to the 
rotational symmetry so it can be easily shown that the 
induced dipole on H2 is of the form in Eq. ([1]) by pro- 
jecting the field perpendicular and parallel to the H2 
bond and calculating the corresponding dipole compo- 
nents. A second contribution to the total dipole of the 
system arises from the H2 permanent quadrupole induc- 



ing a dipole on the MOF. The quadrupolar potential and 
the corresponding electric field at position r, depend on 
r, the H2 quadrupole and the bond orientation, which 
are again of the form in Eq. (pj. This field shifts the 
MOF charge density and induces a dipole. The total 
dipole on the MOF may be formally calculated by multi- 
plying the electric field by the polarizability at the same 
position and integrating over the whole MOF. This pro- 
cedure extends the classical picture of point charge into 
the continuous charge density regime. It cannot be per- 
formed in practice since the polarizability is not available. 
However, the final result for the dipole would be like the 
expression in Eq. (p} , since the integration runs over the 
MOF space while (9, 0) would be left unchanged. One 
can similarly add second-order corrections where the in- 
duced dipole on H 2 and MOF further produce dipole on 
each other. The final equation after this correction turns 
out to be quite simple for cup site absorption (see the 
Supplemental Material for derivation) and reads 

u s x = C{ sin 26> cos - (Cf cos 20 - CJ sin 20) sin 2 6 
My = C{ sin 26> sin + (C^ sin 20 + C3 cos 20) sin 2 9 (2) 
u s z = C% cos 2 9 + CI, 

where s could be H2, MOF, or the total system. The C 
coefficients depend on the H2 quadrupole, polarizability 
and MOF geometry which are kept fixed during the vi- 
brational transition. From Eq. © we see that only two 
orientations (each orientation gives three equations) are 
required to determine the five constants and correspond- 
ingly the dipole for any other orientations. To test this 
model, we calculate u' z for several H2 orientations. Good 
linearity is obtained between u' z and cos 2 9 (see Fig. S7 
in Supplemental Material), in agreement with our model. 
X and Y components are also consistent with our model 
(see Supplemental Material). 

Table Ull summarizes our results for H 2 at the cup site. 
First we consider pure vibrations where rotational quan- 
tum numbers do not change (the Q lines in Table [TTJ) . 
The angular integral for Q(0) (para) is much smaller 
than that for Q(l) (ortho), owing partly to the vanishing 
of the X and Y components of the dipole for the Q(0) 
transition due to symmetry. This symmetry issue also 
applies to the \jm >= |10)— >|10) transition in Q(l) so 
that the integral is about 1 /3 of that for the other two 
transitions of Q(l). Additionally, Fig. [2] shows that para 
H2 has a larger probability to be oriented in XY plane, 
giving a smaller u' z upon bond stretching, while the 1 10) 
state of ortho H2 is p z like and the H 2 bond is mainly 
perpendicular to XY plane. As a result, m'J for the para 
state is only about one quarter of that for the 1 10) state 
(sec Supplemental Material for the integral results of each 
component of u' for the Q transitions). To get the rela- 
tive intensity between Q(0) and Q(l), we need to consider 
the population ratio between para and ortho hydrogen, 
which we took to be 1:3. Also, the calculated rotational 
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TABLE II: Theory vs. experiment [f| for RV transitions at 
the cup site. The frequency shift Av (cm^ 1 ) is relative to 
the corresponding free H2 value. The theoretical intensity 
[oc 1^ times the 30K Boltzmann factor (times 3 for ortho)] is 
normalized to 100 for the strongest line; strong (str), weak 
(wk), and absent (abs) describe the experimental intensity. 





in, 


mf 


Theory 


Experiment 




Av 


Int. 


Av 


Int. 


Q(o) 








-23 


2 




abs 


Q(l) 


±1 



±1 



-23 


97 


-27.5 


str 


Q*(l) 


±1 





22 


9 


39 


wk 






±2 


-44 


58 


-49.3 


str 


S(0) 





±1 


-12 


5 


-6.8 


wk 









-1 


2 




abs 






±3 


-34 


100 


-36.8 


str 




±1 


±2 


-9 


6 


-0.8 


wk 




±1 


6 


9 


21.6 


wk 


S(l) 







11 


3 




abs 




±3 


-78 







abs 







±2 


-53 


3 


-61 


wk 




±1 


-50 


~0 




abs 









-33 


~0 




abs 



energy of the 1 10) state is about 5.5 meV higher than that 
for the m = ±1 states (see Supplemental Material). As 
such, its population is about 13% of that of the m = 1 or 
— 1 state at the experimental temperature of 30K within 
a Boltzmann distribution approximation. We can there- 
fore estimate that the vibrational intensity for para H2 
is about 2.5% of that for ortho H2, hence agreeing with 
the IR measurement, where the para line was simply not 
observed 

TablelTllalso shows the results for S lines in the IR spec- 
trum, where Aj — 2. First a selection rule of Am = ±2 
is observed with small probabilities for other transitions. 
This table also predicts that there is one single strong 
line for each S(0) and S(l) at the experimental tempera- 
ture 30K, with shifts of —44 and —34 cm^ 1 respectively, 
whereas the S(l) line at —53 cm -1 should be weak due to 
the low population of the 1 10) state. More importantly, 
the strong line in each category exhibits the largest fre- 
quency relative to the free H2 value. These results are 
in very good agreement with the IR measurements in 
Ref.| where a single strong S(0) line of —49.3 cm 1 and 
a strong S(l) line of —36.8 cm -1 are observed for H2 at 
the cup site. A weak S(0) line at —6.8 cm -1 and two 
S(l) peaks at —0.8 and 21.6 cm -1 are also observed with 
intensities roughly one order of magnitude smaller than 
that of the corresponding strong line, consistent with our 
calculations. Furthermore, Table HI] shows that the cal- 
culated intensities of the two strong S lines and the Q(l) 
lines are comparable, which is also observed (20| . We 
also note a peak of ^—61 cm -1 shift with an intensity 
similar to that of the S(l) at -0.8 cm' 1 [20]. This peak 
might arise from the |10)— ^ |3, ±2) transitions with theo- 



retical intensity close to those of the other two weak S(l) 
lines, after the 13% population weight is taken into ac- 
cout (Table ITT|). Finally, we discuss the special Q*(l) line 
(|1,±1) — > 1 10) ) ) that is experimentally observed 0. 
The calculated l\ of this transition is approximately 
equal to that for Q(0). However, the population between 
para and ortho hydrogen is ~1:3 which makes Q*(l) 3~4 
times stronger and observable. The calculated shift of 22 
cm -1 is quite small compared to the experimental value 
of 39 cm -1 . This is likely due to the neglect of rotation- 
translation coupling, which would probably lower the 
low rotational state even more and therefore increase the 
splitting between the m=0 and m=±l states. 

In summary, we have proposed a method that pro- 
vides an intuitive picture of H2 interaction in complex 
environments. These techniques provide powerful tools 
for studying gas adsorption in general. 

Supported by DOE Grant No. DE-FG02-08ER46491. 
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Calculational methods 

First-principles calculations based on van der Waals density functional theory were performed within the plane- 
wave implementation of the density functional theory in the ABINIT package which we have adapted from the 
Siesta code to incorporate the van der Waals interaction. We adopted Troullier-Martins pseudopotentials [ij with 
a gradient-corrected functional. An energy cutoff of 50 Ry and Gamma point sampling were used for total energy 
calculations. 



Vibrational frequency 

The four type of adsorption sites are shown in Fig. IS 1 1 To calculate the stretch frequency for H2 at each of these 
four sites, we performed a series of calculations varying the bond length of H2, with the center of H2 and the host 
atoms fixed at their equilibrium positions. The resulting potential-energy curve was used in the Schrodinger equation 
to obtain the eigenvalues and excitation frequencies. A similar calculation was also carried out for isolated H2 to 
obtain the frequency shift due to MOF-H2 interaction. The ab initio total energies vs H2 internuclear distance were 
tabulated in Tables [S3l — IS6I for H2 at the four type of adsorption sites. 



Rotational frequency 

In order to calculate the rotational states, we first sample the solid angle to get the total energies. The spherical 
surface was sampled as follows: the polar angle were evenly divided into seven layers; and the azimuthal angle were 
then sampled by {1,8,16,24,16,8,1} number of points corresponding to each layer from pole to pole of the sphere. We 
next fit these potential energies with spherical harmonics 

V{e,<t>)=Y,ClmYlm (SI) 
im 

which was then substituted into the rigid rotor equation and diagonalized for rotational energies. We found that 
fitting with s and d states gave results converged within 1 cm -1 . The fitted coefficients are shown in Table [S7l The 
calculated rotational energy states are shown in Tables IS8HS111 



Wannier function approach for dipole moment and IR intensity 

The Wannier functions were calculated with the Wannier90 code Q embedded in ABINIT and the Brillouin zone 
was sampled by a 2x2x2 Monkhorst-Pack grid. The Wannier centers for the bare MOF were first calculated and used 
as initial guess for the H2 loaded system. The change of MOF Wannier centers upon H2 adsorption was obtained 
from 

5r n = rf ° F + H * - rf OF (S2) 

where r„ is the center of the n-th Wannier function. Fig. IS3H S5lshow these changes for the 02, 03 and benzene sites 
while the cup site is given in the main text. 

To calculate the IR intensity, one first needs to get the derivative of the dipole with respect to the normal coordinates 
corresponding to H2 stretch vibration. We used the H2 internuclear distance as an approximation for the stretching 
normal coordinates and the derivative was approximate by the finite difference. To reduce the numerical errors, the 
bond stretching should be sufficiently large, but still in the linear regime. We found that a stretch of 0.05A from 
equilibrium bond length was appropriate, as shown in Fig. IS6I 
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Model for induced dipole 

The total induced dipole of the system can be approximated by a sum of four terms. 

u = < 2 + u MOi? + uf OF + uf 2 (S3) 

The first term on the right-hand side is the induced dipole on H 2 due to interactions with MOF atoms. The second 
term is the induced dipole on MOF atoms due to H 2 quadrupole. The third term is the induced dipole on MOF 
due to u^ 2 and the fourth term is the induced dipole on H 2 due to uff OF . These last two terms are second order 
corrections. We now derive the expressions for the four terms for cup site adsorption. 

u^ 2 

Due to the 3-fold symmetry, the electric field (E) at cup site due to MOF atoms is along the rotation axis, i.e. Z. 
For H 2 with its bond oriented along (9, </>), the projected fields along and perpendicular to the bond are 

E\\=E cos 9 (S4a) 

Ej_=Esm6 (S4b) 

and the corresponding induced dipole is 

U|l = Ea\\ cos9 (S5a) 

u± = Eaj_ sin 9 (S5b) 

where a\\ and a± are the H 2 polarizability along and perpendicular to the bond. In Cartesian coordinates, this gives 

— E( a \\ — a ±) sin 9 cos 9 cos (ft (S6a) 
u^y — E(a\\ — a±) sin 9 cos 9 sin (f> (S6b) 
u^J = E{a\\ cos 2 9 + a ± sin 2 9) (S6c) 

,MOF 
U() 

Hydrogen molecule has permanent quadrupole. The tensor is Q zz — —2Q XX = —2Q yy = Q. For H 2 at origin with 
orientation of (9,(j)), the quadrupole potential at position Pi(X,Y,Z) is 

(ST) 

where r = (X 2 + Y 2 + Z 2 ) 1 / 2 and Z = X sin 9 cos <f> + Y sin 9 sin <j> + Z cos 9. The electric field of this potential is 

E QX (r) = |^ {-2r 2 sin9 cos (t>Z + hXZ 2 } - ^X (S8a) 
E 0Y (r) = |^ {-2r 2 sin0sin^ + 5rZ 2 } - (S8b) 
E ox (r) = ^ {-2r 2 sin9 cos <j>Z + 5ZZ 2 } - ^Z (S8c) 

The MOF charge density will be shifted by this field and thus leads to induced dipole. To calculate this induced 
dipole on MOF, one may view that there is an electron at position Pi with partial charge which is equal to the charge 
density at Pi. This charge has a certain polarizability. The total induced dipole can then be formally calculated 
by multiplying the electric held by the corresponding polarizability and then integrating over the whole MOF space. 
This procedure is somewhat an extension of the classical point charge into the continuous charge density regime. 
Assuming the polarizability is isotropic, the final result will have the same dependence on (9, <f)) as electric field since 
the integration runs over (X,Y,Z) while (6, 4>) will be left unchanged. In other words, we will have an equation of the 
form in Eq. (5) in the main text. Note that the isotropic assumption is not critical here except in making the final 
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equations simpler. If one had used the whole polarizability tensor, the final result can still be cast into the form in 
Eq. (5) in the main text. We found that the isotropic assumption gave consistent results for our system. 

Due to the 3-fold rotation symmetry, there are three equivalent points with equal polarizability in MOF. Taking 
advantage of this symmetry, the sum of the electric fields at the three positions is 



E QX = ^ | (3r 2 - hZ 2 ) Z sin 29 cos - ^ (3XY 2 - X 3 ) sin 2 9 cos 20 + | (3X 2 Y - Y 3 ) sin 2 9 sin 20 J 



(S9a) 

E 0Y = ^{ (3r 2 - 5Z 2 ) Z sin 2(9 sin + ^ (3X 2 F - Y 3 ) sin 2 6» cos 20 + ^ (3XY 2 - X 3 ) sin 2 6 sin 20j> (S9b) 

fe.g g^lg pcrf,-!) (S9c) 

The induced dipole on MOF is therefore given by 

U^OF = QMOF gin 2 ^ CQS ^ _ c MOF sin 2 fl cog + gMOF ^2 g gin ^ 

ti M y Of = C ^f OF sin 29 sin + C™ OF sin 2 cos 20 + C^ OF sin 2 sin 20 (SlOb) 

u £|°f = C^ OF (3cos 2 9 - 1) (SlOc) 



where 



c MOF = J 90 _ ^ Za MOF (r) dr ( g Ua) 

C MQF = J ^ 5 _ ^ a MOF (r) dr (snb) 

c MOF^ JW5_ {3x 2 Y _ Y3)a MOF {r)dr (SUc) 

y ^(5Z 2 - 3r 2 )Za MOF (r) rfr (Slid) 



°04 



MOF 



and the integration runs over the 1/3 irreducible region of MOF, as a result of the 3-fold rotation symmetry. 



2 nd -order corrections: uf OF and uf 2 

The induced dipole on H2 and MOF, u^ 2 and u^ fOF , further induces dipole on each other and produces second 
order corrections. For u^ 2 , it gives electric field at Pi(X,Y,Z) 



E 1X (r) = ^ [-u^r 2 + 3< 2 • rX) (S12a) 
E 1Y (r) - 1 (-<^r 2 + 3< 2 • r y) (S12b) 
Siz(r) - ^ (-<|r 2 + 3< 2 • vZ] (S12c) 



Similar to the derivation of u^ OF , one obtains 



u ^ OF = C$ OF sin 20 cos - C^ OF sin 2 9 cos 20 + C^ OF sin 2 9 sin 20 (S13a) 
«* OF = C^ OF sin 20 sin + C^ OF sin 2 9 sin 20 + C^ OF sin 2 9 cos 20 (S13b) 
u m° f = C^ OF cos 2 + eg OF (S13c) 
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where 



riMOF 



riMOF 
°12 



9Q f 3Z _ 5Z^ | 3£^1 3Z 2 



MOF 
13 



riMOF 
°14 



MOF 
15 



4 \ 7* u 7* 1 

/ 27g 

1 4 

9Q /3Z 5Z 3 



3Z + 5Z 3 



y MOF (r)dr 
t MOF (r)dr 
-3E 



, MOF 



(r)dr 



1 



3Z 

5 



3£ 



1 



3Z 2 

^5" 



■) (a||-a±)U MOF (r)dr 



, MOF 



(r) dr 



(S14a) 
(S14b) 
(S14c) 
(S14d) 
(S14e) 



and the integration again runs over 1/3 of the MOF region. 

Now let us look at the second-order correction on H2, u^ 2 . The hydrogen quadrupole generates electric field at Pi 
as given by Eq. (|S9|) . With the help of the partial charge and local polarizability concept, this field produces a local 
dipole UQ /OF (r) = a MOF (r)E (r) where E (r) is given in Eq. (1551) and a MOF (r) is assumed to be isotropic. The 
electric field back at H2 due to this local dipole is 



E 



h 2 

IX 



a 

a 



{E^ OF (r 2 - 3X 2 ) - 3E^ OF XY - 3E^ z OF XZ} 



E $ = ~ {E$° F (r 2 - 3Y 2 ) - 3E^ OF XY - 3E$ z OF YZ} 



E H ' 2 



{E^ z OF (r 2 - 3Z 2 ) - 3E^ x OF XZ - 3E™ OF YZ} 



(S15a) 
(S15b) 
(S15c) 



Inserting Eq. (|S9[) into Eq. (|S15I) . adding together the three rotationally equivalent points and integrating over the 
1 /3 MOF region, one finally obtains 



E ix = C ii sin 29 cos § - C i2 sin2 cos 2( t> + C i3 sin2 sin 2< /> 
e iy = C u sin 29 cos <t> + c i2 sin2 6 sin 2 <t> + C\i sm2 cos 20 



r?Fli 



C^ 2 cos 2 



c 



where 



°12 



r*Fl2 

°13 — 

rtFli 

°14 — 



i)a MOF {r)QZ (33_Z?__ 15Z^ _ 5X 2 Y 2 
16 ~ 8r 2 ~ 16r 4 ~ 



2r» 
9a MOF 



4r 4 



dr 



(r)Q 



2r w 

g a MOF 



(3Y 2 



X 2 )Xdr 



(r)Q 



2r io 
9a MOF 



(3X 2 - Y 2 )Ydr 



(r)QZ 



2r s 



9 45 Z 2 15 Z 4 5X 2 F 2 



16 8 r 2 



16 r 4 



4r 4 



dr 



H 2 
5 



9a MOF (r)QZ ^57 37 Z 2 15 Z 4 5X 2 F 2 
"§"72" ~ Y6 7T 4 r 4 



2r 8 



16 



r/r 



(S16a) 
(S16b) 
(S16c) 



(S17a) 
(S17b) 
(S17c) 
(S17d) 
(S17e) 



and the integral is over the 1/3 of the MOF space. Considering the anisotropy of the polarizability of H2, we have 

(S18) 



H 2 



sin 2 9 cos 2 6 sin 2 9 sin 6 cos 6 sin 9 cos 9 cos < 



a_L + (a|| — a±) sin sin </> cos <fi sin #sin 
sin 9 cos # cos 6 sin cos sin < 



sm cos V sm < 
cos 2 



E 



H 2 



The anisotropic term impose a small correction to the first term inside the bracket. For simplicity, we neglect the 

uf 2 and Ef 2 have £ 
dependence on (9, </>) as given in Eq. f|S16[) . 



second term so that u^ 2 and E^ 2 have a simple linear relationship. In particular, they have the same form of 
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Coefficients 



From Eq. ([56]). (|S10|) . (|S13[) and (IS16|) . we conclude that the following equations hold 




C{ sin 29 cos - (C 2 S cos 20 - C 3 S sin 20) sin 2 6 
C{ sin 26* sin + (CJ sin 20 + CJ cos 20) sin 2 
C 4 S cos 2 6> + CI, 



(S19a) 
(S19b) 
(S19c) 



where s denotes the system and could be H2, MOF or the total. To determine the coefficients C's, we calculated the 
dipole and the derivative of the dipole with respect to H2 internuclear distance with Wannier function approach for 
five hydrogen orientations. The Z components of the obtained values were used to fit C4 and C5 in Eq. (IS19c[) . As 
shown in Fig. IS7[ good linearity is obtained in agreement with our model. 

To compute Ci, C2 and C3, we pick three ab initio calculated value, VL x /u' y of orientation 4 and of orientation 
5, to solve a 3x3 linear equation for the coefficients. To check the values obtained, we substitute them back into 
Eq. (ISf 9al) and (IS19bl) for other orientations and compare with the ab initio results. The comparison are shown in 
Table IS12I and IS131 Consistent results are obtained generally while we do see some deviations on the induced dipole 
on H2, which may be due to the neglect of the anisotropy in Eq. (|S18[) . However, the absolute magnitude of these 
deviations are quite small (< 10%) compared to the total value which is the sum of the induced dipole on MOF and 
on H2. 
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FIG. SI: Illustration of H2 adsorption sites in MOF-5 unit cell (top) and primitive cell (bottom). MOF-5 has FCC structure 
with lattice constant of 25.89A jfj]. The primitive cell has 106 atoms. The H atoms on benzene rings are omitted for clarity. 
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FIG. S2: Illustration of the MOF frame of reference. Origin is at the cup adsorption sites and Z is along the <lll>direction 
of cubic crystal lattice shown in Fig. IS1I It is also the 3-fold rotation axis. 




FIG. S3: Illustration of the 02 adsorption site and the change of Wannier centers due to H2 adsorption compared to bare MOF 
and free H2. The vector lengths are enlarged by 1200. 




FIG. S4: Illustration of the 03 adsorption site and the change of Wannier centers due to H2 adsorption compared to bare MOF 
and free H2. The vector lengths are enlarged by 1200. 




FIG. S5: Illustration of the benzene adsorption site and the change of Wannier centers due to H2 adsorption compared to bare 
MOF and free H2. The vector lengths are enlarged by 1200. 



H 2 




0.02 0.04 0.06 0.08 0.1 



AR(A) 
MOF 




AR(A) 

FIG. S6: Su as a function of H2 internuclear distance, uo is the dipole at equilibrium distance. 



TABLE S3: Calculated total energy vs H2 internuclear distance at cup site 
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TABLE S4: Calculated total energy vs H2 internuclear distance at 02 site 
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TABLE S5: Calculated total energy vs H2 internuclear distance at 03 site 
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2.33191 
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TABLE S6: Calculated total energy vs H2 internuclear distance at benzene site 
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TABLE S7: Expansion coefficients (meV) of orientational potential energy surface in the basis of spherical harmonics. The 
equilibrium energy is set to be zero. 



site 


s 


d z 2 




dyz 


dxy 


d x 2 —y 2 


cup 


17.0 


0.06 


8.45 


8.41 


8.55 


-0.006 


02 


25.0 


-3.74 


-6.36 


-6.90 


-5.76 


-4.88 


03 


7.73 


0.27 


-2.43 


-2.39 


-2.38 





benzene 


1.46 


-0.52 





0.27 





0.80 



TABLE S8: Rotational eigen energies (meV) at cup site 



State # 


Enerffv 


E -Ei 


i 


m 


1 


4.428 




o 





2 


17.498 


13.070 




_1 


3 


17.555 


13.127 


1 


1 


4 


23.014 


18.586 




o 


5 


46.197 


41.769 




_2 


g 


46.197 


41.769 




2 


7 


50.094 


45.666 


2 


_1 


8 


50.135 


45.707 




1 


9 


51.781 


47.353 




o 


10 


89.88 


85.452 




-3 


11 


89.88 


85.452 




3 


12 


92.934 


88.506 




-2 


13 


92.934 


88.506 


3 


2 


14 


94.856 


90.428 




-1 


15 


94.894 


90.466 




1 


16 


95.552 


91.124 








TABLE S9: Rotational eigen energies (meV) at 02 site 



State # 


Enernrv 


E-Ei 


1 


6.76 




2 


18.614 


11.854 


3 


22.306 


15.546 


4 


24.03 


17.270 


5 


49.036 


42.276 


6 


49.39 


42.630 


7 


50.619 


43.859 


g 


53 965 


46.505 


9 


53.439 


46.679 


10 


93.026 


86.266 


11 


93.146 


86.386 


12 


94.315 


87.555 


13 


95.211 


88.451 


14 


95.488 


88.728 


15 


97.788 


91.028 


16 


97.801 


91.041 



TABLE S10: Rotational eigen energies (meV) at 03 site 



State # 


Eneruv 


E -Ei 


i 


m 


1 


2.148 




o 


o 


2 


15.815 


13.667 




o 


3 


17.356 


15.208 


1 


-1 


4 


17.434 


15.286 




1 


5 
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43.403 




o 


g 
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_1 


7 
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43.776 


2 


1 


8 
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44.878 
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2 
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87.736 




1 


13 
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88.227 


3 


-2 


14 


90.377 


88.229 




2 


15 


91.253 


89.105 




-3 


16 


91.253 


89.105 




3 



rpt nr n C11 

lAhSLli oil: 


Rotational eigen energ 


pes (meV) at benzene site 


State # 


Energy 


Ej-Ei 


1 


0.409 


- 


2 


14.929 


14.520 


3 


15.051 


14.642 


4 


15.35 


14.941 


5 


44.332 


43.923 


6 


44.339 


43.930 


7 


44.553 


44.144 


8 


44.64 


44.231 


9 


44.691 


44.282 


10 


88.408 


87.999 


11 


88.409 


88.000 


12 


88.595 


88.186 


13 


88.611 


88.202 


14 


88.693 


88.284 


15 


88.773 


88.364 


16 


88.787 


88.378 
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FIG. S7: duz/dR vs cos 2 9 for the induced dipole on H2 and MOF. u' z was calculated as the difference between the dipole 
moments at equilibrium bond length and a stretch of 0.05A. 
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TABLE S12: Comparison of (8u x , 8u Y ) on MOF due to H 2 bond stretching of 0.05A between the ab initio and the fitted values 



orientation # (degree) 4> (degree) 



ab initio fitted 



5u x 8u y 8u x 8u y 



1 99.650 -95.723 5.690E-04 9.440E-04 5.690E-04 9.634E-04 

2 87.054 -6.190 -4.864E-04 -8.400E-04 -5.178E-04 -8.720E-04 

3 10.097 -112.811 -5.760E-05 -1.008E-04 -5.172E-05 -9.052E-05 

4 54.741 -98.278 2.774E-04 2.464E-04 

5 42.937 -16.295 2.640E-04 -6.126E-04 -5.754E-04 



TABLE S13: comparison of (8ux, 8uy) on Hb due to bond stretching of 0.05A between the ab initio and the fitted values 



orientation 


# (degree) 


^(degree) 


ab initio 


fitted 




8u x 


SUy 


8u x 


SUy 


1 


99.650 


-95.723 


-1.631E-05 


-4.322E-05 


-2.286E-05 


-6.552E-05 


2 


87.054 


-6.190 


1.562E-05 


2.220E-05 


1.340E-05 


4.522E-05 


3 


10.097 


-112.811 


1.338E-05 


3.052E-05 


9.392E-06 


2.030E-05 


4 


54.741 


-98.278 


-2.190E-06 


3.340E-05 






5 


42.937 


-16.295 


-6.386E-05 


3.326E-05 




4.184E-05 
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TABLE S14: Theory vs experiment (Ref.[§j) for RV frequency shifts (cm 1 ) of adsorbed H2 relative to free H2. The vibrational 
transition is from the ground state to the 1st excited state (v=0 — > v=l). 



site 




S(0) (para) 
C?=0-> 3=2) 


S(l) (ortho) 


02 


Th. 


-37, -34, -24, -3, -1 


-15, -14, -4,3, 5, 24, 24 


Ex. 


-36.7, -27.3, -24.3, -7.4 


-12.9 



03 Th. -19, -16, -7 -10, -9, -4, 3 



TABLE S15: Angular integral (jm\u'\jm) for H2 at cup site. The energy of the |00) is set as reference. The units for u' is 
10" 3 e. 



jm 




U'y 




u'^xurV) 


E(meV) 


00 


-0.009 


0.02 


-2.23 


5.0 





11 


-2.38 


7.69 


-1.41 


66.7 


13.1 


11 


2.36 


-7.67 


-1.40 


66.4 


13.2 


10 


0.003 


0.01 


-4.67 


21.8 


18.6 
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TABLE S16: Theoretical predictions and experimental data [f| for v = — > v = 1 transitions of H2 at the cup site. The 
frequency shift Av (cm -1 ) is relative to the corresponding free H2 value and the angular integral given by l\ = \{j /m /|u' \jimi)\ 2 
(10~ 6 e 2 ). The rotational energy (meV) E\ ot of the 1 00) state is set as a reference. The theoretical intensity is calculated from 
l\ weighted by the 30K Boltzmann factor and the spin ratio of 1:3 between para and ortho H2. The strongest line is normalized 
to 100. 





111, 
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Theory 




Experiment 
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Av 


T' A 
1 A 


Intensity 


Av 


Intensity 


Q(0)tf 1 =0-». j f =0) 
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5 


2 




abspnt 


Q(i)0i=i^j/=i) 


±1 




±1 




13.1 
18.6 


-23 


66 
22 


97 


-27.5 


strong 


Q*(l)0'i=l-M/=1) 


±1 





13.1 


22 


6 


9 


39 


weak 






±2 




-44 


115 


58 


-49.3 


strong 


S(0) u. o •./■ 2) 





±1 


13.1 


-12 


10 


5 


-6.8 


weak 











-1 


5 


2 




absent 






±3 




-34 


69 


100 


-36.8 


strong 




±1 


±2 


13.1 


-9 


4 


6 


-0.8 


weak 




±1 


6 


6 


9 


21.6 


weak 


S(l)(ji=l-»J>=3) 









11 


2 


3 




absent 




±3 




-78 










absent 







±2 


18.6 


-53 


49 


3 


-61 


weak 




±1 


-50 


8 


~0 




absent 











-33 


6 


~0 




absent 
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